#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

#define do_divergence_free FALSE
#define ytop_constant (2.0)
#define use_blasius FALSE
#define spread_divide (0.3)

#define DD (1.0)
#define stokes_tolerance (1e-8)
#define fraction_wind (1.05)
#define new_fraction_wind (0.0)
/* above was 0.30, nondim c=.76884=dimensional 3.8m/s */
/* u^2=lg  l=height g=9.8m/s^2 so for above example
   with l=2.5meters,g=9.8 u=4.94 so c_dimensional=u(c_nondim)=3.8.
   for breaking wave, l=2pi g=9.8=>u=7.85 so c_dim=u(c_non)=
     (7.85)(0.43)=3.4    */
/* to get 2m/s wind, we have fraction=2/3.8=0.53 */
/* to get 4m/s wind, we have fraction=4/3.8=1.05 */
/* to get 5m/s wind, we have fraction=5/3.8=1.32 */
/* to get 6m/s wind, we have 6/3.8=1.58 */
/* to get 7m/s wind, we have 7/3.8=1.84 */
/* to get 8m/s wind, we have 8/3.8=2.11 */
/* to get 9m/s wind, we have 9/3.8=2.37 */
/* to get 10m/s wind, we have 10/3.8=2.63 */
/* to get 0m/s wind, we have 0/3.8=0.00 */
/* to get 6.0m/s wind going AGAINST the wave, we have -6/3.8=-1.58 */
/* more notes on dimensionality:
if reference length is L=2.5 meters length of domain is 10 meters=4L, height
of full domain is 5 meters=2L.  Mean height is L.
reference velocity is U^2=Lg  g=9.8m/s^2 L=2.5m
so U=sqrt(2.5*9.8)=4.95
so, we have T=L/U=2.5/4.95=0.51
so, we have t_nondim=t_dim/T
if the step size is 0.05 in non-dim units, then
stepsize is 0.05*T=0.05*0.51=0.0255seconds


For one meter wave we have: l=0.25m g=9.8 u=1.57 T=L/U=0.159

*/
#define dim_wind (5.0)
/* dim_wind in m/s */
/* the fraction_wind stuff has wind travelling in SAME direction as the
   wave */

extern void init_stream_function();

static double nondim_height;
static int stokes_mx;
static double stokes_period;
static double zvalue[2*imax+4];
static int random_k[40];
static double wind_fraction;

double scosh(xx,yy)
double xx,yy;
{
double zz;

  zz=(1.0-exp(-2.0*xx))/(exp(yy-xx)+exp(-yy-xx));
  return zz;
}

double ccosh(xx,yy)
double xx,yy;
{
double zz;

  zz=(1.0+exp(-2.0*xx))/(exp(yy-xx)+exp(-yy-xx));
  return zz;
}

double um_val_gen(xval,yval,zvalue,wavenum)
double xval,yval;
double zvalue[],wavenum;
{
int ix;
double total;

  total=zvalue[stokes_mx+1];
  for (ix=1;ix<=stokes_mx;ix++)
   total+=wavenum*ix*zvalue[stokes_mx+1+ix]*
     ccosh(ix*wavenum*yval,ix*wavenum*DD)*cos(ix*pi*xval);
  return total;
}

double um_val(ival,zvalue,wavenum)
int ival;
double zvalue[],wavenum;
{
double tempx;

tempx=ival;
tempx/=stokes_mx;
return um_val_gen(tempx,zvalue[ival],zvalue,wavenum);
}

double vm_val_gen(xval,yval,zvalue,wavenum)
double xval,yval;
double zvalue[],wavenum;
{
int ix;
double total;

  total=0.0;
  for (ix=1;ix<=stokes_mx;ix++)
   total+=wavenum*ix*zvalue[stokes_mx+1+ix]*
     scosh(ix*wavenum*yval,ix*wavenum*DD)*sin(ix*pi*xval);
  return total;
}

double vm_val(ival,zvalue,wavenum)
int ival;
double zvalue[],wavenum;
{
double tempx;

tempx=ival;
tempx/=stokes_mx;
return vm_val_gen(tempx,zvalue[ival],zvalue,wavenum);
}

double sjm_val(index,jval,mval,zvalue,wavenum)
double zvalue[],wavenum;
int index,jval,mval;
{
double factor;

  factor=scosh(jval*wavenum*zvalue[mval],jval*wavenum*DD);
  if (index==1)
   return factor*cos(jval*mval*pi/stokes_mx);
  else
   return factor*sin(jval*mval*pi/stokes_mx);
}


double cjm_val(index,jval,mval,zvalue,wavenum)
double zvalue[],wavenum;
int index,jval,mval;
{
double factor;

  factor=ccosh(jval*wavenum*zvalue[mval],jval*wavenum*DD);
  if (index==1)
   return factor*sin(jval*mval*pi/stokes_mx);
  else
   return factor*cos(jval*mval*pi/stokes_mx);
}

void get_fvalue(zvalue,fvalue,height,wavenum,residual)
double zvalue[],fvalue[],height,wavenum,*residual;
{
int ix,iy;
double um_temp,vm_temp;

  for (ix=0;ix<=stokes_mx;ix++)
   {
    fvalue[ix]=zvalue[ix]*zvalue[stokes_mx+1]+zvalue[2*stokes_mx+2];
    for (iy=1;iy<=stokes_mx;iy++)
     fvalue[ix]+=zvalue[stokes_mx+1+iy]*
      scosh(iy*wavenum*zvalue[ix],iy*wavenum*DD)*cos(iy*ix*pi/stokes_mx);

    um_temp=um_val(ix,zvalue,wavenum);
    vm_temp=vm_val(ix,zvalue,wavenum);
    fvalue[ix+stokes_mx+1]=(sqr(um_temp)+sqr(vm_temp))/2.0+zvalue[ix]-
     zvalue[2*stokes_mx+3];
   }

  fvalue[2*stokes_mx+2]=(zvalue[0]+zvalue[stokes_mx])/(2.0*stokes_mx)-1;
  for (ix=1;ix<=stokes_mx-1;ix++)
   fvalue[2*stokes_mx+2]+=zvalue[ix]/stokes_mx;

  fvalue[2*stokes_mx+3]=zvalue[0]-zvalue[stokes_mx]-height;

  *residual=0.0;
  for (ix=0;ix<=2*stokes_mx+3;ix++)
   if (fabs(fvalue[ix])>*residual)
    *residual=fabs(fvalue[ix]);
}   /* end get_fvalue */

double get_stokes_height(zvalue,ix)
double zvalue[];
int ix;
{
int ix_alt,iz;
double xval,total,stokes_h,xj;

  ix_alt=ix;
  if (2*ix>mx) ix_alt=mx-ix;

  stokes_h=h*mx/(stokes_mx*2.0);

  xval=ix_alt*h;   /* NO xmin here, xmin=0 in paper */
  total=0.0;
  for (iz=0;iz<=stokes_mx;iz++)
   {
    xj=iz*stokes_h;
    if ((xval>xj-stokes_h)&&(xval<=xj))
     total+=zvalue[iz]*(1.0+(xval-xj)/stokes_h);
    else if ((xval>=xj)&&(xval<=xj+stokes_h))
     total+=zvalue[iz]*(1.0-(xval-xj)/stokes_h);
   }
  return total;
}   /* get_stokes_height */

double elevation(ix)
int ix;
{
int iy;
double hold,f1,f2;

 hold=(ymax+ymin)/2.0;
 for (iy=0;iy<=my;iy++)
  {
   f1=d1[ix][iy];f2=d1[ix][iy+1];
   if (((f1<0.0)&&(f2>=0.0))||
       ((f1>=0.0)&&(f2<0.0)))
    hold=ymin+(iy-0.5)*h+h*fabs(f1)/(fabs(f1)+fabs(f2));
  }

 return hold;
}   /* elevation */

void speed_interface(ix,iy,uinterface,vinterface,zvalue,wavenum)
int ix,iy;
double *uinterface,*vinterface;
double zvalue[],wavenum;
{
int ix_alt;
double elevate,xval,yval;

 ix_alt=ix;
 if (2*ix>mx) ix_alt=mx-ix;
 xval=h*(ix_alt-0.5);   /* NO xmin here, xmin=0 in paper */
 xval*=(wavenum/pi);  /* stupid fix here */
 yval=ymin+h*(iy-0.5);
 elevate=elevation(ix);
 if (yval>elevate)
  yval=elevate;

 *uinterface=um_val_gen(xval,yval,zvalue,wavenum);
 *vinterface=vm_val_gen(xval,yval,zvalue,wavenum);
 if (2*ix>mx) *vinterface=-(*vinterface);
}   /* speed_interface */
  
void u_from_s()
{
int ix,iy;

 for_domain
  {
   u[ix][iy]=(s[ix][iy]-s[ix][iy-1]+s[ix-1][iy]-s[ix-1][iy-1])/(2.0*h);
   v[ix][iy]=-(s[ix][iy]-s[ix-1][iy]+s[ix][iy-1]-s[ix-1][iy-1])/(2.0*h);
  }
}
  
/* checks d1 for current height, updates u,v,s,top_stream_value  */
void stokes_add_wind()
{
int ix,iy;
double change_wind,yval,end_point,height_min,height_max,height_value,
 velocity_spread,stream_max,stream_min,stream_value,yroof,yfloor;
int iyfloor;

 for_all_points
  if (! point_ok)
   {
    d1[ix][iy]=gddd(d1,ix,iy);
    s[ix][iy]=gst1(s,ix,iy);
   }

 velocity_spread=nondim_height/spread_divide;
 change_wind=(fraction_wind-new_fraction_wind)*zvalue[stokes_mx+1];
 /* above is (1-new)-(1-old) */

 for (ix=0;ix<=mx;ix++)
 {
  height_value=(elevation(ix)+elevation(ix+1))/2.0;
  iyfloor=(int) ((height_value-ymin)/h);
  yfloor=ymin+iyfloor*h;
  yroof=yfloor+h;
  stream_value=(yroof-height_value)*s[ix][iyfloor]/h +
    (height_value-yfloor)*s[ix][iyfloor+1]/h;

  if (ix==0)
   { stream_min=stream_max=stream_value;
     height_min=height_max=height_value; }
  else
   {
    if (stream_min<stream_value) stream_min=stream_value;
    if (stream_max>stream_value) stream_max=stream_value;
    if (height_min<height_value) height_min=height_value;
    if (height_max>height_value) height_max=height_value;
   }
 
  for (iy=0;iy<=my;iy++)
   {
    yval=ymin+iy*h;
    if (yval>height_value)
     {
      s[ix][iy]+=(change_wind/2.0)*(yval-height_value+
        velocity_spread*log(
          cosh((yval-height_value)/velocity_spread ) ) );

      if (new_fraction_wind==0.0)
       s[ix][iy]=stream_value+zvalue[stokes_mx+1]*(yval-height_value);
     }
   }  /* looping iy */
 }  /* looping ix */
 top_stream_value+=(change_wind/2.0)*(ymax+fary-1.0+
   velocity_spread*log(
     cosh((ymax+fary-1.0)/velocity_spread ) ) );
 if (new_fraction_wind==0.0)
  {
   printf("stream_min,max height_min,max %10.5f %10.5f %10.5f %10.5f\n",
     stream_min,stream_max,height_min,height_max);
   top_stream_value=(stream_min+stream_max)/2.0+
      zvalue[stokes_mx+1]*(ymax+fary-(height_min+height_max)/2.0);
  }

 printf("new top_stream_value %10.5f \n",top_stream_value);

 u_from_s();
}   /* stokes_add_wind */

double get_stokes_stream(zvalue,wavenum,ix,iy)
double zvalue[],wavenum;
int ix,iy;
{
int ix_alt,iz;
double xval,yval,height_value,total,start_point,end_point,
 velocity_spread,bigval,yval1;
double nondim_speed,Qval,ycritical,kappa,alpha,yparam,yzero,
 utop_fraction,ytemp,coeff3,coeff4,uint;

  ix_alt=ix;
  if (2*ix>mx) ix_alt=mx-ix;

  xval=ix_alt*h;   /* NO xmin here, xmin=0 in paper */
  yval=ymin+iy*h;
 
  height_value=get_stokes_height(zvalue,ix);
  velocity_spread=nondim_height/spread_divide;
  if (stokes_breaking_flag)
   velocity_spread/=40.0;

  if (stokes_wind_flag) wind_fraction=wind_parameter;
  nondim_speed=fabs(zvalue[stokes_mx+1]);
  Qval=zvalue[2*stokes_mx+2];
  ycritical=1.0+1.25*nondim_height;
  ycritical=1.0;
  kappa=0.41;
  alpha=0.0144;
  yparam=alpha*sqr(wind_fraction*nondim_speed);

  if ((ix==0)&&(iy==0)&&(stokes_wind_flag)&&
      (use_blasius))
   {
    yzero=ycritical+yparam*
     (exp(kappa/wind_fraction)-1.0);
    utop_fraction=(wind_fraction/kappa)*log(1.0+
     (2.0-ycritical)/yparam)-1.0;
    printf("|c|,ycrt,wdfrc,yzero,utpfrc%10.5f%10.5f%10.5f%10.5f%10.5f\n",
     nondim_speed,ycritical,wind_fraction,yzero,utop_fraction);
   }

#define fake_height 1.0E6
#define wind_done_below FALSE

  if ((yval<=height_value)||
      ((wind_done_below)&&(stokes_wind_flag)) )
  {
   total=-nondim_speed*yval;
   if (yval<=height_value)
    yval1=yval;
   else if ((yval<=height_value+nondim_height+h)&&(stokes_wind_flag))
    yval1=yval;
   else
    yval1=height_value*(fake_height-yval)/(fake_height-height_value);

   for (iz=1;iz<=stokes_mx;iz++)
    total+=zvalue[stokes_mx+iz+1]*
     scosh(iz*wavenum*yval1,iz*wavenum*DD)*cos(iz*wavenum*xval);
  }
  else
  {
   if (yval<=ytop_constant)
    {
     uint=um_val_gen((xval*wavenum/pi),height_value,zvalue,wavenum);
     ytemp=ytop_constant-height_value;
     coeff4=(uint-nondim_speed)/sqr(ytemp)+
       2.0*nondim_speed*(ytop_constant-1.0)/(sqr(ytemp)*ytemp);
     coeff3=(-nondim_speed-uint-3.0*coeff4*sqr(ytemp))/(2.0*ytemp);
     ytemp=yval-height_value;
     total=-Qval+uint*ytemp+coeff3*sqr(ytemp)+coeff4*sqr(ytemp)*ytemp;
    }
   if (((stokes_wind_flag)||(stokes_breaking_flag))&&(wind_fraction!=0.0))
    {
     total=-Qval-nondim_speed*(yval-1.0);

     if (yval>ycritical)
      {
       if (use_blasius)
        {
         ytemp=1.0+(yval-ycritical)/yparam;
         total+=(nondim_speed*wind_fraction*yparam/kappa)*
            (ytemp*(log(ytemp)-1.0)+1.0);
        }
       else
        {
         bigval=(yval-1.0)/velocity_spread;
         if (bigval>=500.0)
          bigval=(yval-1.0)-velocity_spread*log(2.0);
         else
          bigval=velocity_spread*log(cosh(bigval));

         total+=wind_fraction*nondim_speed*bigval;
        }
      }   /* start wind above mean level */

    }   /* doing wind case */
     
  }  /* updating air */

  return total;
}  /* get_stokes_stream */

void matrix_solve(AA,xx,bb)
double AA[2*imax+4][2*imax+4],xx[],bb[];
{
double AA_save[2*imax+4][2*imax+4],bb_save[2*imax+4],
  alpha,hold_value,hold_sum,bb_sum,AA_sum,hold1;
int ix,iy,iz,NN,hold_iy;

  NN=2*stokes_mx+3;
bb_sum=0.0;
for (ix=0;ix<=NN;ix++)
 bb_sum+=fabs(bb[ix]);
printf("bb sum %e \n",bb_sum);

AA_sum=0.0;
  for (ix=0;ix<=NN;ix++)
   {
    bb_save[ix]=bb[ix];
    for (iy=0;iy<=NN;iy++)
{
AA_sum+=fabs(AA[ix][iy]);
     AA_save[ix][iy]=AA[ix][iy];
}
   }
printf("AA sum %e \n",AA_sum);
  
  for (ix=0;ix<=NN-1;ix++)
   {
    hold_iy=ix;hold_value=AA[ix][ix];
    for (iy=ix+1;iy<=NN;iy++)
     if (fabs(AA[iy][ix])>hold_value)
      {
       hold_iy=iy;
       hold_value=fabs(AA[iy][ix]);
      }

    if (hold_iy!=ix)
     for (iy=ix;iy<=NN;iy++)
      {
       hold_value=AA[ix][iy];
       AA[ix][iy]=AA[hold_iy][iy];
       AA[hold_iy][iy]=hold_value;
      }
    hold_value=bb[ix];
    bb[ix]=bb[hold_iy];
    bb[hold_iy]=hold_value;

    for (iy=ix+1;iy<=NN;iy++)
     {
      alpha=AA[iy][ix]/AA[ix][ix];
      for (iz=ix;iz<=NN;iz++)
       AA[iy][iz]-=alpha*AA[ix][iz];
      bb[iy]-=alpha*bb[ix];
     }
   }   /* gaussian elimination w/partial pivoting */ 

  for (ix=NN;ix>=0;ix--)
   {
    hold_value=bb[ix];
    for (iy=ix+1;iy<=NN;iy++)
     hold_value-=AA[ix][iy]*xx[iy];
    xx[ix]=hold_value/AA[ix][ix];
   }

  hold_value=0.0;
  hold_sum=0.0;
  for (ix=0;ix<=NN;ix++)
   {
    hold_sum+=fabs(xx[ix]);
    hold1=bb_save[ix];
    for (iy=0;iy<=NN;iy++)
     hold1-=AA_save[ix][iy]*xx[iy];
    hold_value+=fabs(hold1);
   }
  printf("matrix solve validation resid,xx sum %e %e\n",
    hold_value,hold_sum);
}    /* matrix_solve */

#define max_wave_numbers 20

double random_profile(xx)
double xx;
{
int ix;
double temp_wavelen,sum_series;

  temp_wavelen=(xmax-xmin);
  sum_series=0.0;
  for (ix=0;ix<max_wave_numbers;ix++)
   sum_series+=(nondim_height/random_k[ix])*
     sin(2.0*pi*random_k[ix]*(xx-xmin)/temp_wavelen);
  return (sum_series+(ymin+ymax)/2.0);
}

/* input: mx,my,dimensionless wavelen and height,xmin,xmax,ymin,ymax,
     fary
   output: s,u,v,d1 

NOTE: we double the mesh, but since psi and eta are symmetric about
 x, we cut in half the number of points we solve for.
*/
init_stokes_data(dim_wavelen,dim_height)
double dim_wavelen,dim_height;
{
double height,wavelen,wavenum,wave_slope,wave_period;
double jacobian[2*imax+4][2*imax+4],fvalue[2*imax+4],
  new_zvalue[2*imax+4],wave_speed,residual,
  c_sum,s_sum,current_wave_slope,dmiddle,ymiddle,nondim_wind;
double stream_a,stream_b,stream_c,stream_crest;
int ix,iy,iz,mx_all,count,stokes_repeat,stokes_again,test_k;
double velocity_spread,xval,yval,bigval,ycritical,uint,vint,
 heavyfn,uwind,kappa,alpha,nondim_speed,yzero,utop_fraction,
 Qval,ytemp,yparam;
logical found;

  if (stokes_random)
   {
    ix=0;
/* for SUN computers, iy=srandom(1) is ok, for gnu compiler, "srandom(1)" is 
   ok.  */
/*    iy=srandom(1);*/  /* seed the random number generator */
/*    srandom(1);   */
    while (ix<max_wave_numbers)
     {
      test_k=(random()&255);
      if ((test_k>=10)&&(test_k<=40))
       {
        found=FALSE;
        for (iy=0;iy<ix;iy++)
         if (random_k[iy]==test_k)
          found=TRUE;
        if (! found)
         {
          random_k[ix]=test_k;
          printf("seeded ix,k %5d %5d \n",ix,test_k);
          ix++;
         }
       }  /* test_k in the range */
     }  /* looping until 20 wavenumbers found */

    wavelen=(dim_wavelen/(2.0*dim_height));  /* narrow box here 10/2.5 */
    nondim_height=wavelen/(30.0*log(max_wave_numbers));
    nondim_height=wavelen/1000.0;
    velocity_spread=wavelen/60.0;
    nondim_wind=wind_speed/sqrt(9.8*2.0*dim_height);

    for_domain
     d1[ix][iy]=y[iy]-random_profile(x[ix]); 
    for_all_points
     if (!point_ok)
      d1[ix][iy]=gddd(d1,ix,iy);

    for (ix=0;ix<=mx;ix++)
     for (iy=0;iy<=my;iy++)
      {
       xval=xmin+ix*h;
       yval=ymin+iy*h;
       s[ix][iy]=0.0;
       ycritical=random_profile(xval);
/*       ycritical=(ymax+ymin)/2.0;*/
       if (yval>ycritical)
        s[ix][iy]=(nondim_wind/2.0)*(yval-ycritical+
           velocity_spread*log(
             cosh((yval-ycritical)/velocity_spread ) ) );
      }
    u_from_s();
    ycritical=(ymax+ymin)/2.0;
    top_stream_value=(nondim_wind/2.0)*(ymax+fary-ycritical+
           velocity_spread*log(
             cosh((ymax+fary-ycritical)/velocity_spread ) ) );
    printf("top_stream_value %10.5f\n",top_stream_value);
    return;
   }   /* if stokes_random */
    

  if (standing_wave_flag)
   {
    wavelen=(dim_wavelen/dim_height);  /* 2 for our case */

    if (stokes_wind_flag)
     {
      wavelen=2.0; 
      /* non-dimensionalize using dim_wavelen/2 not dim_height*/

      /* was dim_wind, now we use wind_speed */
      nondim_wind=wind_speed/sqrt(9.8*dim_wavelen/2.0);
     }

    wavenum=(2.0*pi)/wavelen;
    wave_slope=0.025;
    if (stokes_wind_flag)
     wave_slope=0.1;

    height=wave_slope/wavenum;
    wave_period=2.0*pi/sqrt(wavenum*tanh(wavenum));
    if (stokes_wind_flag)
     wave_period=2.0*pi/sqrt(wavenum);  /* infinite height */
    wave_speed=wavelen/wave_period;

    for_domain
     d1[ix][iy]=y[iy]-height*cos(wavenum*x[ix]);  /* no more +pi */
    for_all_points
     if (!point_ok)
      d1[ix][iy]=gddd(d1,ix,iy);
/*
    stream_b=nondim_wind/2.0;
    stream_a=stream_b/(2.0*height);
    stream_c=-stream_a*sqr(height)+stream_b*height;
    stream_crest=stream_a*sqr(height)+stream_b*height+stream_c;
*/
    for (ix=0;ix<=mx;ix++)
     for (iy=0;iy<=my;iy++)
      {
       s[ix][iy]=0.0;
       ymiddle=(y[iy]+y[iy+1])/2.0;
       dmiddle=(d1[ix][iy]+d1[ix+1][iy]+d1[ix][iy+1]+d1[ix+1][iy+1])/4.0;
       if ((stokes_wind_flag)&&(dmiddle>=0.0))
        {
         s[ix][iy]=nondim_wind*dmiddle;
/*
         if (ymiddle<=height)
          s[ix][iy]=stream_a*sqr(ymiddle)+stream_b*ymiddle+stream_c;
         else
          s[ix][iy]=stream_crest+nondim_wind*(ymiddle-height);
*/
        }
      }    
     
    u_from_s();  

    top_stream_value=nondim_wind*(ymax+fary);

    printf("standing wave ht,period,k %10.5f %10.5f %10.5f \n",
     height,wave_period,wavenum);
    printf("wave_speed,wave_slope,nondim_wind %10.5f %10.5f %10.5f \n",
     wave_speed,wave_slope,nondim_wind);
    printf("top_stream_value %10.5f\n",top_stream_value);
    return;
   }
     
    
  if ( (mx/2)*2 != mx) printf(" mx must be divisible by 2 \n");

  wavelen=(dim_wavelen/dim_height);  /* dimensionless form */
  wavenum=(2.0*pi)/wavelen;
  wave_slope=0.1;   /* already dimensionless */
  wind_fraction=fraction_wind;
  if (stokes_wind_flag) wind_fraction=wind_parameter;

  stokes_mx=32;   /* was mx */
  if (stokes_breaking_flag)
   {
    wind_fraction=wind_speed;
    wave_slope=0.4;
   }
  else
   wave_slope=wind_speed;
  printf("wind fraction is %10.5f \n",wind_fraction);

  if (wave_slope<=0.0)
   stokes_mx=mx;    /* need the accuracy here */

printf("WAVE SLOPE IS THE WIND SPEED FOR NON BREAKING %10.5f\n",wave_slope);


#define max_wave_slope 0.25
#define diff_wave_slope 0.02

  if (wave_slope<max_wave_slope)
   stokes_again=1;
  else
   stokes_again=(int) ((wave_slope-max_wave_slope)/diff_wave_slope)+2;

for (stokes_repeat=0;stokes_repeat<stokes_again;stokes_repeat++)
{
if (stokes_repeat==stokes_again-1)
 current_wave_slope=wave_slope;
else 
 current_wave_slope=max_wave_slope+stokes_repeat*diff_wave_slope;

  if (current_wave_slope>wave_slope)
   current_wave_slope=wave_slope;

  height=2.0*current_wave_slope/wavenum;  /* peak to peak height */
  nondim_height=height;                 /* peak to peak */
 
  wave_speed=sqrt(tanh(wavenum)/wavenum);
  wave_period=2*pi/(wavenum*wave_speed);


if (stokes_repeat==0)
{
  mx_all=2*stokes_mx+3;
  for (ix=0;ix<=stokes_mx;ix++)
   {
    zvalue[ix]=1.0+height*cos(ix*pi/stokes_mx)/2.0;
    zvalue[ix+stokes_mx+1]=0.0;
   }
  
  zvalue[stokes_mx+1]=-wave_speed;
/*  zvalue[stokes_mx+2]=-height/(4.0*wave_speed*wavenum);*/
  zvalue[stokes_mx+2]=(height/(2.0*wave_speed*wavenum))*ccosh(DD,1.0);

  zvalue[2*stokes_mx+2]=wave_speed;  /* initial value for Q */
  zvalue[2*stokes_mx+3]=1.0+sqr(wave_speed)/2.0;  /* initial value for P */
}

  get_fvalue(zvalue,fvalue,height,wavenum,&residual);

  printf("initial,ht,k,wavelen,c,resid %10.5f%10.5f%10.5f%10.5f%10.5f\n",
   height,wavenum,wavelen,wave_speed,residual);
  printf("wave period 2*pi/(ck) %10.5f\n",wave_period);
  printf("stokes_repeat,current_wave_slope %5d %10.5f \n",stokes_repeat,
   current_wave_slope);

  count=1;

  while ((residual>stokes_tolerance)&&(! stokes_linear_flag))
   {
    
    for (ix=0;ix<=mx_all;ix++)
     {
      for (iy=0;iy<=mx_all;iy++)
       jacobian[ix][iy]=0.0;
      if (ix<=stokes_mx)
       {
        jacobian[ix][ix]=um_val(ix,zvalue,wavenum);
        jacobian[ix][stokes_mx+1]=zvalue[ix];
        for (iy=stokes_mx+2;iy<=2*stokes_mx+1;iy++)
         jacobian[ix][iy]=sjm_val(1,iy-stokes_mx-1,ix,zvalue,wavenum);
        jacobian[ix][2*stokes_mx+2]=1.0;
       }
      else if (ix<=2*stokes_mx+1)
       {
        s_sum=c_sum=0.0;
        for (iz=1;iz<=stokes_mx;iz++)
         {
          s_sum+=sqr(iz)*
            zvalue[stokes_mx+iz+1]*sjm_val(1,iz,ix-stokes_mx-1,zvalue,wavenum);
          c_sum+=sqr(iz)*
            zvalue[stokes_mx+iz+1]*cjm_val(1,iz,ix-stokes_mx-1,zvalue,wavenum);
         }
        jacobian[ix][ix-stokes_mx-1]=1.0+
         um_val(ix-stokes_mx-1,zvalue,wavenum)*sqr(wavenum)*s_sum+
         vm_val(ix-stokes_mx-1,zvalue,wavenum)*sqr(wavenum)*c_sum;
        jacobian[ix][stokes_mx+1]=um_val(ix-stokes_mx-1,zvalue,wavenum);
        for (iy=stokes_mx+2;iy<=2*stokes_mx+1;iy++)
         jacobian[ix][iy]=
           (iy-stokes_mx-1)*wavenum*um_val(ix-stokes_mx-1,zvalue,wavenum)*
           cjm_val(2,iy-stokes_mx-1,ix-stokes_mx-1,zvalue,wavenum)+
           (iy-stokes_mx-1)*wavenum*vm_val(ix-stokes_mx-1,zvalue,wavenum)*
           sjm_val(2,iy-stokes_mx-1,ix-stokes_mx-1,zvalue,wavenum);
        jacobian[ix][2*stokes_mx+3]=-1.0;
       }
      else if (ix==2*stokes_mx+2)
       {
        jacobian[ix][0]=jacobian[ix][stokes_mx]=1.0/(2.0*stokes_mx);
        for (iy=1;iy<=stokes_mx-1;iy++)
         jacobian[ix][iy]=1.0/stokes_mx;
       }
      else if (ix==2*stokes_mx+3)
       {
        jacobian[ix][0]=1.0;
        jacobian[ix][stokes_mx]=-1.0;
       }
     }    /*  end looping ix */

    matrix_solve(jacobian,new_zvalue,fvalue);  /* computes z=Ainv*(f) */
    for (ix=0;ix<=mx_all;ix++)
     zvalue[ix]-=new_zvalue[ix];  /* minus sign because we solve z=-Ainv(f) */
    get_fvalue(zvalue,fvalue,height,wavenum,&residual);
    count++;

    printf("count,resid %5d %20.10f\n",count,residual);
   }  /* while residual not met */
}  /* looping successive heights */

nondim_speed=fabs(zvalue[stokes_mx+1]);
Qval=zvalue[2*stokes_mx+2];
velocity_spread=nondim_height/spread_divide;
if (stokes_breaking_flag)
 velocity_spread/=40.0; 
if (use_blasius)
 {
  ycritical=1.0;
  kappa=0.41;
  alpha=0.0144;
  yparam=alpha*sqr(wind_fraction*nondim_speed);
 }
 
top_stream_value=-Qval-nondim_speed*(ymax+fary-1.0);

if ((wind_fraction!=0.0)&&
    ((stokes_wind_flag)||(stokes_breaking_flag)))
 {
  bigval=(ymax+fary-1.0)/velocity_spread;
  printf("bigval %10.5e\n",bigval);
  if (bigval>=500.0)
   bigval=(ymax+fary-1.0)-velocity_spread*log(2.0);
  else
   bigval=velocity_spread*log(cosh(bigval));
  printf("updated multiplier %10.5e\n",bigval);
  if (use_blasius)
   {
    ytemp=1.0+(ymax+fary-ycritical)/yparam;
    bigval=(yparam/kappa)*(ytemp*(log(ytemp)-1.0)+1.0);
   }

  top_stream_value+=wind_fraction*nondim_speed*bigval;
 }

bottom_stream_value=0.0;     
stokes_period=wavelen/nondim_speed;

 for_all_points
  d1[ix][iy]=y[iy]-(get_stokes_height(zvalue,ix)+
          get_stokes_height(zvalue,ix-1))/2.0;

 for (ix=0;ix<=mx;ix++)
  for (iy=0;iy<=my;iy++)
   s[ix][iy]=get_stokes_stream(zvalue,wavenum,ix,iy);
 u_from_s();   /* stream function version always available */

 if (! do_divergence_free)
  {   
   /* first initialize the steady state solution */
   for_domain
    {
     speed_interface(ix,iy,&uint,&vint,zvalue,wavenum);
     yval=y[iy];
     if (yval>ytop_constant)
      {
       vint=0.0;
       uint=zvalue[stokes_mx+1];
      }
     else if (d1[ix][iy]>0.0)
      {
       heavyfn=sinh(ytop_constant-yval)/sinh(ytop_constant-yval+d1[ix][iy]);
       vint*=heavyfn;
       uint=heavyfn*uint-(1.0-heavyfn)*(nondim_speed);
      }

     if (wind_fraction != 0.0)  /* do not override for non-forced case */
      {
       u[ix][iy]=uint;
       v[ix][iy]=vint;
      }
    }   /* initialize steady solution */

   /* add the wind if non-zero */
   if (wind_fraction!=0.0)
    {
     for_domain
      {
       yval=y[iy];
       if (d1[ix][iy]>0.0)
        {
         heavyfn=tanh((yval-1.0)/velocity_spread);
         if (use_blasius)
          heavyfn=(1.0/kappa)*log(1.0+(yval-ycritical)/yparam );

         if (yval>=1.0)
          u[ix][iy]=-nondim_speed+(wind_fraction)*nondim_speed*heavyfn;
         else
          u[ix][iy]=-nondim_speed;
         v[ix][iy]=0.0;
        }   /* updating air speed */
      }  /* looping all points */
    }   /* adding wind */

if (FALSE)
{
   init_stream_function();  /* uses top_stream and bottom_stream ... */
   for_all_points
    if (iy<iymin)
     s[ix][iy]=bottom_stream_value;  /* should be zero */
    else if (iy>iymax)
     s[ix][iy]=top_stream_value;   /* some big negative number */
    else if (ix<ixmin)
     s[ix][iy]=s[ixmax][iy];
    else if (ix>ixmax)
     s[ix][iy]=s[ixmin][iy];
    else
     s[ix][iy]=pcg_result[ix][iy];
   u_from_s();
}
  }   /* if not divfree */

printf("B0,Q,top_stream %10.5f  %10.5f  %10.5f\n",
  zvalue[stokes_mx+1],zvalue[2*stokes_mx+2],top_stream_value);
printf("period (c=-B0),(lambda/c) %10.5f \n",stokes_period);

}   /* init_stokes_data */

double compute_surface_force(ix)
int ix;
{
double factor,xx,tt,hold;

  factor=2.0*pi/(xmax-xmin);
  if (ix<=0)
   ix=ix+mx;
  else if (ix>mx)
   ix=ix-mx;

  xx=(x[ix]-xmin)*factor;
  tt=2.0*pi*rk_time/stokes_period;
/* surface tension coeff already divided by factor */
  if (tt<=pi)
   hold=sin(tt)*sin(xx-tt+pi);
  else
   hold=0.0;
/*
  if (hold != 0.0)
  printf("hold,ix,xx,tt %10.5f %5d %10.5f %10.5f \n",hold,ix,xx,tt);
*/
  return hold;
}
